!     ---------------
!     PROGRAMME JANUS   GF-9812
!     ---------------
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des coordonnees de la Lune en utilisant les donnees de:
!     Lunar Tables Programs from 4000 B.C. to A.D. 8000
!     M. Chapront-Touze et J. Chapront (Willmann Bell, 1991).
!
!     Les tables de la Lune sont lues sur le fichier : janus.dat.
!     Les resultats sont ecrits sur le fichier impression : janus.txt.
!     Ces deux fichiers sont dans le meme repertoire que le programme.
!
!V98  REMARQUE:
!V98  LA VERSION 1998 DIFFERE DE LA VERSION INITIALE.
!V98  ON NE DISTINGUE PLUS LES 3 NIVEAUX DE PRECISION.
!V98  ON PREND SYSTEMATIQUEMENT LA PRECISION LA PLUS HAUTE (IPR=1).
!V98  LES MODIFICATIONS SONT SIGNALEES PAR 'V98'. 
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character rep
      character*30 datcal
!
! --- Initialisation ---------------------------------------------------
!
      data d2000/2451545.d0/
      data sjul/36525.d0/
!
      nultab=4
      nulprn=10
!
      ipr=1
      icd=1
      ir=1
      idt=1
      idf=2
      d=2451545.0d0
      datcal='01 Janv 2000   12h00m00s'
      ndat=10
      pas=1.0d0
!
! --- Fichier Impression -----------------------------------------------
!
      open (nulprn,file='janusfr.txt')
      write (nulprn,1000)
!
! --- Lecture des Tables de la Lune ------------------------------------
!
      call LECTAB (1,1,nultab)
!
! --- Lecture des Parametres -------------------------------------------
!
100   continue
!
      call MENU (ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
!V98  call LECTAB (2,ipr,nultab)
!
! --- Debut du calcul --------------------------------------------------
!
      t0=d
      dt=0.d0
      pasj=pas/24.d0
!
      select case (icd)
!
! --- Calcul des coordonnees de la Lune --------------------------------
!
      case (1,2)
!
      do i=1,ndat
!
         if (mod(i,10).eq.1) call TITRE (1,nulprn,ipr,icd,ir,idt,idf)
         if (i.eq.1) call TITRE (2,nulprn,ipr,icd,ir,idt,idf)
!
         dj=t0+(i-1)*pasj
         t=(dj-d2000)/sjul
!
         if (idt.eq.2) then
            call UT (dj,dt,ierr,limut,tlimut)
            t=t+dt/315576.d4
         endif
!
         call LUNPOL (t,c1,c2,c3,ir,icd,ipr)
!
         select case (ir)
         case (1)
            call IMPC (dj,nulprn,c1,c2,c3,ipr,dt,idt,idf)
         case (2)
            call IMPQ (dj,nulprn,c1,c2,c3,ipr,dt,idt,idf)
         case (3)
            call IMPC (dj,nulprn,c1,c2,c3,ipr,dt,idt,idf)
         case (4)
            call IMPQ (dj,nulprn,c1,c2,c3,ipr,dt,idt,idf)
         endselect
!
         if (mod(i,10).eq.0) then
            write (*,'(/2x,a,$)') 'Appuyez sur Entree'
            read (*,*)
         endif
!
      enddo
!
! --- Calcul des elements orbitaux de la Lune --------------------------
!
      case (3)
!
      do i=1,ndat
!
         if (mod(i,10).eq.1) call TITRE (1,nulprn,ipr,icd,ir,idt,idf)
         if (i.eq.1) call TITRE (2,nulprn,ipr,icd,ir,idt,idf)
!
         dj=t0+(i-1)*pasj
         t=(dj-d2000)/sjul
!
         call LUNOSC (t,a,e,gam,pom,gom,alam)
         call IMPS (dj,nulprn,a,e,gam,pom,gom,alam,ipr,idf)
!
         if (mod(i,10).eq.0) then
            write (*,'(/2x,a,$)') 'Appuyez sur Entree'
            read (*,*)
         endif
!
      enddo
!
! --- Fin du calcul ----------------------------------------------------
!
      endselect
!
! --- Un autre calcul ? ------------------------------------------------
!
      write (*,1001)
      read (*,1002) rep
!
      if (rep.eq.'n'.or.rep.eq.'N') then
         write (nulprn,'(//)')
         stop
      endif
!
      goto 100
!
! --- Formats ----------------------------------------------------------
!
1000  format (2x,'JANUS GF9812',17x,'LUNE de -4000 a 8000',17x,
     .        'M&J CHAPRONT'/2x,78('-'))
1001  format (/4x,'Le calcul est termine, ',
     .        'voulez-vous effectuer un autre calcul (o/n) ? ',$)
1002  format (a)
!
      end
!
!
!
      BLOCK DATA
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
!     Termes seculaires de la Longitude moyenne de la Lune.
!
! ----------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      common/long/al(5)
!
      data al/
     .218.31665d0,481267.88134d0,-13.268d-4,1.856d-6,-1.534d-8/
!
      end
!
!
!
      subroutine MENU (ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Affichage du Menu.
!
! --- Entree -----------------------------------------------------------
!
!     ipr      Indice Precision (entier).
!     icd      Indice Coordonnees (entier).
!     ir       Indice Repere (entier).
!     idt      Indice Temps (entier).
!     idf      Indice Date Julienne / Date Calendrier (entier).
!     d        Date Julienne initiale (double reel).
!     datcal   Date Calendrier initiale (character*30).
!     ndat     Nombre de dates (entier).
!     pas      Pas du calcul en heures (double reel).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character rep
      character*15 z
      character*30 apr(3),acd(3),arp(4),adt(2),datcal
!
! --- Initialisation ---------------------------------------------------
!
      data apr/
     .'Haute precision',
     .'Moyenne precision',
     .'Faible precision'/
      data acd/
     .'Coordonnees Geometriques',
     .'Coordonnees Apparentes',
     .'Elements orbitaux osculateurs'/
      data arp/
     .'Repere Ecliptique Moyen',
     .'Repere Equatorial Moyen',
     .'Repere Ecliptique Vrai ',
     .'Repere Equatorial Vrai '/
      data adt/
     .'Temps Dynamique (ET)',
     .'Temps Universel (UT)'/
!
! --- Afficher le Menu -------------------------------------------------
!
100   continue
!
      call EFFECR
      write (*,1000)
!
!V98  write (*,1001) apr(ipr)
      write (*,1002) acd(icd)
      write (*,1003) arp(ir)
      write (*,1004) adt(idt)
!
      if (idf.eq.1) write (*,1005) d
      if (idf.eq.2) write (*,1006) datcal
!
      write (z,'(i15)') ndat
      i=1
      do while (z(i:i).eq.' ')
         i=i+1
      enddo
      write (*,1007) z(i:15)
!
      write (z,'(f15.4)') pas
      i=1
      do while (z(i:i).eq.' ')
         i=i+1
      enddo
      if (z(i:i).eq.'.') then
         i=i-1
         z(i:i)='0'
      endif
      write (*,1008) z(i:15)
!
      ierr=0
      write (*,1009)
      read (*,'(a)') rep
!
      select case (rep)
!V98  case ('P','p')
!V98     call LECPRM (1,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('C','c')
         call LECPRM (2,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('R','r')
         if (icd.eq.3) then
            ir=1
            write (*,2001)
            read (*,*)
         else
            call LECPRM (3,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
         endif
      case ('T','t')
         if (icd.eq.3) then
            idt=1
            write (*,2002)
            read (*,*)
         else
            call LECPRM (4,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
         endif
      case ('D','d')
         call LECPRM (5,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('N','n')
         call LECPRM (6,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('I','i')
         call LECPRM (7,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('X','x')
         return
      case ('Q','q')
         stop
      case default
         ierr=1
      endselect
!
      if (ierr.eq.1) then
         write (*,2003)
         read (*,*)
      endif
!
      goto 100
!
! --- Formats ----------------------------------------------------------
!
1000  format (2x,'JANUS GF9812',17x,'LUNE de -4000 a 8000',17x,
     .        'M&J CHAPRONT'/2x,78('-')//
     .        8x,'LISTE DES COMMANDES',t40,'PARAMETRES DU CALCUL'/)
!V98  1001 format (8x,'P : Precision     ',t40,a)
1002  format (8x,'C : Coordonnees    ',t40,a)
1003  format (8x,'R : Repere         ',t40,a)
1004  format (8x,'T : Temps          ',t40,a)
1005  format (8x,'D : Date Initiale  ',t40,'DJ',f13.5)
1006  format (8x,'D : Date Initiale  ',t40,a)
1007  format (8x,'N : Nombre de dates',t40,a)
1008  format (8x,'I : Intervalle (heure)',t40,a)
1009  format (8x,'X : Executer le calcul'/
     .        8x,'Q : Quitter le programme'//
     .        8x,'Entrez votre choix ? ',$)
!
2001  format (/8x,'Pour les Elements orbitaux : ',
     .        'Repere Ecliptique Moyen de la date'/
     .        /8x,'Appuyez sur Entree ',$)
2002  format (/8x,'Pour les Elements orbitaux : ',
     .        'Temps Dynamique (ET)'/
     .        /8x,'Appuyez sur Entree ',$)
2003  format (/8x,'Erreur : Appuyez sur Entree ',$)
!
      end
!
!
!
      subroutine LECTAB (ilec,ipr,nul)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Lecture des tables de la Lune: Longitude, Latitude et Distance.
!     
! --- Entree -----------------------------------------------------------
!
!     ilec     Indice de lecture (entier).
!              ilec = 1 : lecture du fichier des tables de la Lune.
!              ilec = 2 : lecture des tables par niveau de precision.
!
!     ipr      Indice de la precision (entier).
!              ipr = 1 pour Haute precision.
!              ipr = 2 pour Moyenne precision.
!              ipr = 3 pour Faible precision.
!
!     nul      Numero d'unite logique du fichier Tables (entier).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*80 zlec(1324)
!
      dimension nc(3,3),np(3,3,3)
      dimension a(5),b(2)
!
      common/cv/deg,hr
      common/lu/pbc(218,3),apbc(218,5,3),per(423,3),aper(423,2,3),
     .          kc(3),kp(3,3)
!
! --- Initialisation ---------------------------------------------------
!
      data pi/3.141592653589793d0/
!
      data nc/218,188,154,59,45,40,29,14,18/
!
      data np/
     .244,64,114,154,64,68,25,12,19,
     .3,2,0,16,6,9,5,0,4,1,0,0,6,0,3,1,0,0/
!
      deg=pi/180.d0
      hr=deg*15
      il=0
!
! --- Lecture du fichier des tables de la Lune -------------------------
!
      if (ilec.eq.1) then
         call EFFECR
         write (*,1000)
         write (*,1001)
         open (nul,file='janus.dat',status='old')
         do i=1,1324
            read (nul,2000) zlec(i)
         enddo
         close (nul)
         write (*,1002)
         read (*,*)
      endif
!
! --- Lecture des tables par niveau de precision -----------------------
!
      do i=1,3
!
         kc(i)=nc(i,ipr)
         jmax=kc(i)
!
         do j=1,jmax
!
            il=il+1
            if (i.ne.3) then
               read (zlec(il),2001) c,(a(j1),j1=1,5)
            else
               read (zlec(il),2002) c,(a(j1),j1=1,5)
               a(1)=a(1)+90.d0
            endif
!
            pbc(j,i)=c
!
            do j1=1,2
               apbc(j,j1,i)=a(j1)*deg
            enddo
!
            apbc(j,3,i)=a(3)*deg*1.d-4
            apbc(j,4,i)=a(4)*deg*1.d-6
            apbc(j,5,i)=a(5)*deg*1.d-8
!
         enddo
!
         jsup=nc(i,1)-jmax
!
         do j=1,jsup
            il=il+1
            read (zlec(il),2003)
         enddo
!
         jmax=0
!
         do k=1,3
!
            kp(i,k)=np(i,k,ipr)
            jmin=jmax+1
            jmax=jmax+kp(i,k)
!
            do j=jmin,jmax
!
               il=il+1
               if (i.ne.3) then
                  if (k.eq.1) read (zlec(il),3001) c,(b(j1),j1=1,2)
                  if (k.eq.2) read (zlec(il),4001) c,(b(j1),j1=1,2)
                  if (k.eq.3) then
                     read (zlec(il),5001) c,(b(j1),j1=1,2)
                     c=c*1.d-4
                  endif
                  c=c*1.d-3
               else
                  if (k.eq.1) read (zlec(il),3002) c,(b(j1),j1=1,2)
                  if (k.eq.2) read (zlec(il),4002) c,(b(j1),j1=1,2)
                  if (k.eq.3) then
                     read (zlec(il),5002) c,(b(j1),j1=1,2)
                     c=c*1.d-4
                  endif
                  b(1)=b(1)+90
               endif
!
               per(j,i)=c
!
               do j1=1,2
                  aper(j,j1,i)=b(j1)*deg
               enddo
!
            enddo
!
            jsup=np(i,k,1)-kp(i,k)
!
            do j=1,jsup
               il=il+1
               read (zlec(il),2003)
            enddo
!
         enddo
!
      enddo
!
      return
!
! --- Formats ----------------------------------------------------------
!
1000  format (///
     .        2x,t20,'         PROGRAMME JANUS GF9812       '//
     .        2x,t20,'               REFERENCE:             '//
     .        2x,t20,'       LUNAR TABLES AND PROGRAMS      '/
     .        2x,t20,'      FROM 4000 B.C. TO A.D. 8000     '//
     .        2x,t20,'   M. CHAPRONT-TOUZE AND J. CHAPRONT  '//)
1001  format (2x,t20,'   Lecture des Tables de la Lune ...  '/)
1002  format (2x,t20,' Lecture terminee : Appuyez sur Entree',$)
!
2000  format (a)
2001  format(5x,f11.8,2x,f11.7,2x,f16.7,2(2x,f8.3),2x,f7.3)
2002  format(5x,f11.4,2x,f11.7,2x,f16.7,2(2x,f8.3),2x,f7.3)
2003  format(1x)
!
3001  format(5x,f8.5,2x,f8.4,2x,f13.4)
3002  format(5x,f8.4,2x,f8.4,2x,f13.4)
!
4001  format(5x,f8.5,2x,f7.3,2x,f12.3)
4002  format(5x,f8.4,2x,f7.3,2x,f12.3)
!
5001  format(5x,f5.2,2x,f5.1,2x,f10.1)
5002  format(5x,f4.1,2x,f5.1,2x,f10.1)
!
      end
!
!
!
      subroutine LECPRM (iprm,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Lecture des parametres.
!
! --- Entree -----------------------------------------------------------
!
!     iprm     Indice des parametre (entier).
!              1 : ipr    Precision.
!              2 : icd    Coordonnees.
!              3 : ir     Repere.
!              4 : idt    Temps.
!              5 : idf    Date Julienne / Data Calendrier.
!              5 : d      Date initiale (Julienne).
!              5 : datcal Date initiale (Calendrier).
!              6 : ndat   Nombre de dates.
!              7 : pas    Intervalle en heures.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*16 adat,apas
      character*30 datcal,datx
!
! --- Initialisation ---------------------------------------------------
!
      djmin=260057.0d0
      djmax=4643366.d0
      datmin=-40000101.0d0
      datmax=80010101.0d0
      call UT (d,dt,ierr,limut,tlimut)
!
! --- Selection des parametres -----------------------------------------
!
      select case (iprm)
!
! --- Precision --------------------------------------------------------
!V98
!V98  case (1)
!V98
!V98  write (*,1001)
!V98  read (*,'(i10)',iostat=nerr) ix
!V98  if (nerr.ne.0) goto 100
!V98  if (ix.lt.1.or.ix.gt.3) goto 100
!V98  ipr=ix
!V98!
! --- Coordonnees ------------------------------------------------------
!
      case (2)
!
      write (*,1002)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.3) goto 100
      icd=ix
      if (icd.eq.3) then
         ir=1
         idt=1
      endif
!
! --- Repere -----------------------------------------------------------
!
      case (3)
!
      write (*,1003)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.4) goto 100
      ir=ix
!
! --- Temps ------------------------------------------------------------
!
      case (4)
!
      write (*,1004)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.2) goto 100
!
      df=d+pas*(ndat-1)/24.d0
      if (ix.eq.2.and.(df-tlimut).gt.0.d0) goto 400
      idt=ix
!
! --- Date initiale ----------------------------------------------------
!
      case (5)
!
      write (*,1005)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.2) goto 100
      idf=ix
!
      if (idf.eq.1) then
         write (*,1006)
         read (*,'(f20.8)',iostat=nerr) x
         if (nerr.ne.0) goto 100
         if (x.eq.0.d0) return
         if (x.lt.djmin.or.x.gt.djmax) goto 100
      else
         write (*,1007)
         read (*,'(a)',iostat=nerr) adat
         if (nerr.ne.0) goto 100
         if (adat.eq.' ') return
         n=index(adat,'.')
         if (n.eq.0) then
            read (adat,'(i16)',iostat=nerr) ix
            if (nerr.ne.0) goto 100
            xx=ix
         else
            read (adat,'(f16.6)',iostat=nerr) xx
            if (nerr.ne.0) goto 100
         endif
         if (xx.lt.datmin.or.xx.gt.datmax) goto 100
         call DCDJ (xx,x,ierr,datx)
         if (ierr.eq.1) goto 100
         if (ierr.eq.2) goto 200
      endif
!
      df=x+pas*(ndat-1)/24.d0
      if ((df-djmax).gt.0.d0) goto 300
      if (idt.eq.2.and.(df-tlimut).gt.0.d0) goto 400
      d=x
      call DJDATE (d,3,ij,im,ia,nh,nm,ns)
      date=abs(ia)*1.d4+im*1.d2+ij+nh*1.d-2+nm*1.d-4+ns*1.d-6
      if (ia.lt.0) date=-date
      call DCDJ (date,x,ierr,datcal)
!
! --- Nombre de dates --------------------------------------------------
!
      case (6)
!
      write (*,1008)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.1000) goto 100
!
      df=d+pas*(ix-1)/24.d0
      if ((df-djmax).gt.0.d0) goto 300
      if (idt.eq.2.and.(df-tlimut).gt.0.d0) goto 400
      ndat=ix
!
! --- Intervalles ------------------------------------------------------
!
      case (7)
!
      write (*,1009)
      read (*,'(a)',iostat=nerr) apas
      if (ierr.ne.0) goto 100
      if (apas.eq.' ') goto 100
      n=index(apas,'.')
      if (n.eq.0) then
         read (apas,'(i16)',iostat=nerr) ix
         if (nerr.ne.0) goto 100
         x=ix
      else
         read (apas,'(f16.8)',iostat=nerr) x
         if (nerr.ne.0) goto 100
      endif
      if (x.le.0.d0.or.x.gt.1000.d0) goto 100
!
      df=d+x*(ndat-1)/24.d0
      if ((df-djmax).gt.0.d0) goto 300
      if (idt.eq.2.and.(df-tlimut).gt.0.d0) goto 400
      pas=x
!
! --- Fin de la selection ----------------------------------------------
!
      endselect
      return
!
! --- Messages Erreur --------------------------------------------------
!
100   continue
!
      write (*,2001)
      read (*,*)
      return
!
200   continue
!
      write (*,2002)
      read (*,*)
      return
!
300   continue
!
      write (*,2003)
      read (*,*)
      return
!
400   continue
!
      write (*,2004) limut
      read (*,*)
      return
!
! --- Formats ----------------------------------------------------------
!
!V98  1001  format (/8x,'1: Haute precision'
!V98       .        /8x,'2: Moyenne precision'
!V98       .        /8x,'3: Faible precision',
!V98       .        t40,'Votre choix ? ',$)
1002  format (/8x,'1: Coordonnees polaires geometriques'
     .        /8x,'2: Coordonnees polaires apparentes'
     .        /8x,'3: Elements orbitaux osculateurs',
     .        t50,'Votre choix ? ',$)
1003  format (/8x,'1: Ecliptique et Equinoxe moyens de la date'
     .        /8x,'2: Equateur et Equinoxe moyens de la date'
     .        /8x,'3: Ecliptique moyen et Equinoxe vrai de la date'
     .        /8x,'4: Equateur et Equinoxe vrais de la date',
     .        t60,'Votre choix ? ',$)
1004  format (/8x,'1: Temps Dynamique (ET)'
     .        /8x,'2: Temps Universel (UT)',
     .        t40,'Votre choix ? ',$)
1005  format (/8x,'1: Dates Julienne'
     .        /8x,'2: Date Calendrier',
     .        t40,'Votre choix ? ',$)
1006  format (/8x,'Entrez une Date Julienne (decimal) : ',$)
1007  format (/8x,'Entrez une Date Calendrier ',
     .        '(+/-aaaammjj.hhmmss) : ',$)
1008  format (/8x,'Entrez le Nombre de dates (entier) : ',$)
1009  format (/8x,'Entrez le Pas du calcul en heures (decimal) : ',$)
!
2001  format (/8x,'Erreur : Appuyez sur Entree ',$)
2002  format (/8x,'Erreur : Date invalide dans le calendier Gregorien',
     .        ' : Appuyez sur Entree ',$)
2003  format (/8x,'Erreur : Date au-dela de 8000',
     .        ' : Appuyez sur Entree ',$)
2004  format (/8x,'Erreur : Date posterieure a la limite UT ',i4,
     .        ' : Appuyez sur Entree ',$)
!
      end
!
!
!
      subroutine LUNPOL (t,c1,c2,c3,ir,icd,ipr)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des coordonnees polaires de la Lune.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
!     ir       Indice du Repere (entier).
!              ir = 1 pour le Repere ecliptique moyen de la date
!              ir = 2 pour le Repere equatorial moyen de la date.
!              ir = 3 pour le Repere ecliptique vrai.
!              ir = 4 pour le Repere equatorial vrai.
!
!     icd      Indice des coordonnees (entier).
!              icd = 1 pour les coordonnees geometriques.
!              icd = 2 pour les coordonnees apparentes.
!
!     ipr      Indice de la precision (entier).
!              ipr = 1 pour Haute precision.
!              ipr = 2 pour Moyenne precision.
!              ipr = 3 pour Faible precision.
!
! --- Output -----------------------------------------------------------
!
!     c1       Longitude (deg) ou Ascension Droite (h) (double reel).
!     c2       Latitude ou Declinaison en degres (double reel).
!     c3       Distance Terre-Lune en kilometres (double reel).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
! --- Calcul -----------------------------------------------------------
!
      select case (ir)
      case (1)
         call MECC (t,c1,c2,c3,icd)
      case (2)
         call MEQC (t,c1,c2,c3,icd)
      case (3)
         call TECC (t,c1,c2,c3,icd,ipr)
      case (4)
         call TEQC (t,c1,c2,c3,icd,ipr)
      endselect
!
      return
      end
!
!
!
      subroutine MECC (t,along,alat,r,icd)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des coordonnees ecliptiques moyenne de la Lune.
!     Longitude et latitude de la Lune, Distance Terre-Lune.
!     Repere : Ecliptique et Equinoxe moyens de la date.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
!     icd      Indice des coordonnees (entier).
!              icd = 1 pour les coordonnees geometriques.
!              icd = 2 pour les coordonnees apparentes.
!
! --- Sortie -----------------------------------------------------------
!
!     along    Longitude de la Lune en degres (double reel).
!     alat     Latitude de la Lune en degres (double reel).
!     r        Distance Terre-Lune en kilometres (double reel).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension c(3)
!
      common/cv/deg,hr
      common/long/al(5)
      common/lu/pbc(218,3),apbc(218,5,3),per(423,3),aper(423,2,3),
     .          kc(3),kp(3,3)
!
! --- Calcul -----------------------------------------------------------
!
      do i=1,3
!
         c(i)=0
!
         do j=1,kc(i)
            phi=apbc(j,5,i)
            do j1=4,1,-1
               phi =phi*t+apbc(j,j1,i)
            enddo
            s=sin(phi)
            c(i)=c(i)+s*pbc(j,i)
         enddo
!
         isup=kp(i,1)+kp(i,2)+kp(i,3)
         ct=0
!
         do it=3,1,-1
            inf=isup-kp(i,it)+1
            cc=0
            do j=inf,isup
               phi=aper(j,1,i)+t*aper(j,2,i)
               s=sin(phi)
               cc=cc+s*per(j,i)
            enddo
            ct=cc+t*ct
            isup=inf-1
         enddo
!
         c(i)=c(i)+ct
!
      enddo
!
      phi=al(5)
!
      do i=4,1,-1
         phi=phi*t+al(i)
      enddo
!
      along=c(1)+phi
      alat=c(2)
      r=385000.57d0+c(3)
!
      if (icd.eq.2) then
         phi=(225.d0+477198.9d0*t)*deg
         along=along-0.19524d-3-0.1059d-4*sin(phi)
         r=r+0.708d-1*cos(phi)
         phi=(183.3d0+483202.d0*t)*deg
         alat=alat-0.1754d-4*sin(phi)
      endif
!
      along=mod(along,360.d0)
      if(along.lt.0) along=along+360
!
      return
      end
!
!
!
      subroutine MEQC (t,ra,de,r,icd)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des coordonnees equatoriales moyennes de la Lune.
!     Ascension droire et Declinaison de la Lune, Distance Terre-Lune.
!     Repere : Equateur et Equinoxe moyens de la date.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
!     icd      Indice des coordonnees (entier).
!              icd = 1 pour les coordonnees geometriques.
!              icd = 2 pour les coordonnees apparentes.
!
! --- Sortie -----------------------------------------------------------
!
!     ra       Ascension droite de la Lune en heures (double reel).
!     de       Declinaison de la Lune en degres (double reel).
!     r        Distance Terre-Lune en kilometres (double reel).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      common/cv/deg,hr
!
! --- Calcul -----------------------------------------------------------
!
      call OBLI (t,epsi)
!
      epsi=epsi*deg
      ce=cos(epsi)
      se=sin(epsi)
!
      call MECC (t,alg,alt,r,icd)
!
      alg=alg*deg
      alt=alt*deg
      ca=cos(alg)
      sa=sin(alg)
      cb=cos(alt)
      sb=sin(alt)
!
      xq1=cb*ca
      xe2=sa*cb
      xq2=ce*xe2-se*sb
      xq3=se*xe2+ce*sb
!
      ra=atan2(xq2,xq1)
      de=asin(xq3)
      ra=ra/hr
      de=de/deg
      if (ra.lt.0) ra=ra+24
!
      return
      end
!
!
!
      subroutine TECC (t,along,alat,r,icd,ipr)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des coordonnees ecliptiques vraies de la Lune.
!     Longitude et latitude de la Lune, Distance Terre-Lune.
!     Repere : Ecliptique moyen et Equinoxe vrai de la date.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
!     icd      Indice des coordonnees (entier).
!              icd = 1 pour les coordonnees geometriques.
!              icd = 2 pour les coordonnees apparentes.
!
!     ipr      Indice de la precision (entier).
!              ipr = 1 pour Haute precision.
!              ipr = 2 pour Moyenne precision.
!              ipr = 3 pour Faible Precision.
!
! --- Sortie -----------------------------------------------------------
!
!     along    Longitude de la Lune en degres (double reel).
!     alat     Latitude de la Lune en degres (double reel).
!     r        Distance Terre-Lune en kilometres (double reel).
!
! --- Calcul -----------------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      call MECC (t,along,alat,r,icd)
      call NUTA (t,dpsi,deps,ipr)
!
      along=along+dpsi
      along=mod(along,360.d0)
      if(along.lt.0) along=along+360
!
      return
      end
!
!
!
      subroutine TEQC (t,ra,de,r,icd,ipr)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des coordonnees equatoriales vraies de la Lune.
!     Ascension droite et Declinaison de la Lune, Distance Terre-Lune.
!     Repere : Equateur et Equinoxe vrais.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
!     icd      Indice des coordonnees (entier).
!              icd = 1 pour les coordonnees geometriques.
!              icd = 2 pour les coordonnees apparentes.
!
!     ipr      Indice de la precision (entier).
!              ipr = 1 pour Haute precision.
!              ipr = 2 pour Moyenne precision.
!              ipr = 3 pour Faible Precision.
!
! --- Sortie -----------------------------------------------------------
!
!     ra       Ascension droite de la Lune en heures (double reel).
!     de       Declinaison de la Lune en degres (double reel).
!     r        Distance Terre-Lune en kilometres (double reel).
!
! --- Calcul -----------------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      common/cv/deg,hr
!
      call OBLI (t,epsi)
      call NUTA (t,dpsi,deps,ipr)
!
      epsi=epsi+deps
      epsi=epsi*deg
      ce=cos(epsi)
      se=sin(epsi)
!
      call MECC (t,alg,alt,r,icd)
!
      alg=alg+dpsi
      alg=alg*deg
      alt=alt*deg
!
      ca=cos(alg)
      sa=sin(alg)
      cb=cos(alt)
      sb=sin(alt)
!
      xq1=cb*ca
      xe2=sa*cb
      xq2=ce*xe2-se*sb
      xq3=se*xe2+ce*sb
!
      ra=atan2(xq2,xq1)
      de=asin(xq3)
      ra=ra/hr
      de=de/deg
      if (ra.lt.0) ra=ra+24
!
      return
      end
!
!
!
      subroutine LUNOSC (t,a,e,gam,pom,gom,alam)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des elements orbitaux osculateurs de la Lune.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
! --- Sortie -----------------------------------------------------------
!
!     a        Demi-grand axe en km (double reel).
!     e        Exentricite (double reel).
!     gam      Sinus de la demi inclinaison (double reel).
!     pom      Longitude du perigee en degres (double reel).
!     gom      Longitude du noeud en degres (double reel).
!     alam     Longitude moyenne en degres (double reel).
!
!     Repere : Ecliptique et Equinoxe moyens de la date.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      common/cv/deg,hr
!
      data g/4.0184172d24/
!
! --- Calcul -----------------------------------------------------------
!
      call DERIV (t,v,u,r,vp,up,rp)
!
      v=v*deg
      u=u*deg
!
      csv=cos(v)
      ssv=sin(v)
      csu=cos(u)
      ssu=sin(u)
!
      x=vp*csu*csu
      y=up*ssv-vp*ssu*csu*csv
      z=up*csv+vp*ssu*csu*ssv
!
      w=sqrt(x*x+y*y+z*z)
      wp=w*deg
      den=sqrt(2*w*(w+x))
      gc=z/den
      gs=y/den
      gam2=gc*gc+gs*gs
      gam=sqrt(gam2)
      gom=atan2(gs,gc)
      gom=gom/deg
      if (gom.lt.0) gom=gom+360.d0
!
      rw=r*wp
      rw2=rw*rw
      a=1.d0/(2.d0/r-(rp*rp+rw2)/g)
      c=ssu/sqrt(1-gam2)
      xp=csu*csv-c*gs
      yp=csu*ssv+c*gc
      zp=rw2*r/g-1
      ws=rp*r*rw/g
      ec=xp*zp+yp*ws
      es=yp*zp-xp*ws
      e2=ec*ec+es*es
      e=sqrt(e2)
      pom=atan2(es,ec)
      pom=pom/deg
      if (pom.lt.0) pom=pom+360.d0
!
      raci=sqrt(1-e2)
      c=ws/(raci*(1+raci))
      rsa=r/a
      cepi=ec+rsa*(xp-es*c)
      sepi=es+rsa*(yp+ec*c)
      epi=atan2(sepi,cepi)
      alam=epi-(ec*sepi-es*cepi)
      alam=alam/deg
      alam=mod(alam,360.d0)
      if (alam.lt.0) alam=alam+360.d0
!
      return
      end
!
!
!
      subroutine OBLI (t,epsi)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Calcul de l'Obliquite moyenne.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
! --- Sortie -----------------------------------------------------------
!
!     epsi     Obliquite moyenne en degres (double reel).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      dimension eps(5)
!
      data eps/23.43928d0,-0.13d-1,0.d0,0.555d-6,-0.14d-9/
!
! --- Obliquite moyenne ------------------------------------------------
!
      epsi=eps(5)
!
      do k=4,1,-1
         epsi=epsi*t+eps(k)
      enddo
!
      return
      end
!
!
!
      subroutine NUTA (t,dpsi,deps,ipr)
! ----------------------------------------------------------------------
!
!     Ref : GF-0512
!
! --- Objet ------------------------------------------------------------
!
!     Calcul de la Nutation en longitude et en obliquite.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
!     ipr      Indice de la precision (entier).
!              ipr = 1 pour Haute precision.
!              ipr = 2 pour Moyenne precision.
!              ipr = 3 pour Faible Precision.
!
! --- Sortie -----------------------------------------------------------
!
!     dpsi     Nutation en longitude en degres (double reel).
!     deps     Nutation en obliquite en degres (double reel).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension dps(13),dep(9),ar(13,3),dps1(5),dep1(2),kn(3)
!
      common/cv/deg,hr
!
! --- Initialisation ---------------------------------------------------
!
      data kn/13,2,1/
!
      data dps/
     . -4.77767d-3,-0.36631d-3,-0.6317d-4, 0.5728d-4,
     .   0.3961d-4,  0.1978d-4,-0.1436d-4,-0.1072d-4, -0.836d-5,
     .    0.603d-5,  -0.439d-5,  0.358d-5, 0.342d-5/
!
      data dps1/
     . -0.484d-5,3*0.d0,-0.9d-7/
!
      data dep/
     . 2.55625d-3,0.15933d-3,0.2714d-4,-0.2486d-4,
     . 2*0.d0,0.622d-5,0.556d-5,0.358d-5/
!
      data dep1/
     . 0.25d-6,-0.9d-7/
!
      data ar/
     .    125.0446d0,   200.9329d0,    76.6333d0,    250.0891d0,
     .    357.5291d0,   134.9634d0,   198.4620d0,    311.5888d0,
     .    211.5967d0,
     .    203.4038d0,   259.2630d0,     75.8883d0,   301.6699d0,
     .  -1934.1362d0, 72001.5397d0, 962535.7627d0, -3868.2724d0,
     .  35999.0503d0,477198.8676d0, 108000.5899d0,964469.8989d0,
     .1439734.6303d0, 36002.4894d0,-413335.3554d0, 73935.6758d0,
     . 485336.8951d0,
     .  21.d-4,6.d-4,-27.d-4, 42.d-4,  -2.d-4,90.d-4,5.d-4,-47.d-4,
     .  63.d-4,8.d-4,123.d-4,-15.d-4,-117.d-4/
!
      dpsi=0
      deps=0
!
! --- Nutation en longitude et en obliquite ----------------------------
!
      do i=1,kn(ipr)
!
         phi=ar(i,3)
!
         do j=2,1,-1
            phi=phi*t+ar(i,j)
         enddo
!
         phi=phi*deg
         c=dps(i)
         if (i.eq.1.or.i.eq.5) c=c+t*dps1(i)
         dpsi=dpsi+c*sin(phi)
!
         if (i.ne.5.and.i.ne.6.and.i.le.9) then
            c=dep(i)
            if (i.eq.1.or.i.eq.2) c=c+t*dep1(i)
            deps=deps+c*cos(phi)
         endif
!
      enddo
!
      return
      end
!
!
!
      subroutine DERIV (t,v,u,r,vp,up,rp)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Calcul des coordonnees moyennes de la Lune et de leurs derivees.
!
! --- Entree -----------------------------------------------------------
!
!     t        Temps en siecles juliens depuis J2000 (double reel).
!
! --- Sortie -----------------------------------------------------------
!
!     v        Longitude en degres (double reel).
!     u        Latitude en degres (double reel).
!     r        Distance en km (double reel).
!     vp       Longitude derivee en degres/siecle (double reel).
!     up       Latitude derivee en degres/siecle (double reel).
!     rp       Distance derivee en km/siecle (double reel).
!
!     Repere : Ecliptique et Equinoxe de la date.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension c(3),dc(3),dal(4)
!
      common/lu/pbc(218,3),apbc(218,5,3),per(423,3),aper(423,2,3),
     .          kc(3),kp(3,3)
      common/long/al(5)
!
      data dal/481267.881d0,-0.26536d-2,0.5568d-5,-0.6136d-7/
!
! --- Calcul -----------------------------------------------------------
!
      do i=1,3
!
         c(i)=0
         dc(i)=0
!
         do j=1,kc(i)
            phi=apbc(j,5,i)
            do j1=4,1,-1
               phi=phi*t+apbc(j,j1,i)
            enddo
            dphi=4*apbc(j,5,i)
            do j1=4,2,-1
               dphi=dphi*t+(j1-1)*apbc(j,j1,i)
            enddo
            c(i)=c(i)+pbc(j,i)*sin(phi)
            dc(i)=dc(i)+pbc(j,i)*dphi*cos(phi)
         enddo
!
         isup=kp(i,1)+kp(i,2)+kp(i,3)
         dct=0
         ct=0
!
         do it=3,1,-1
            inf=isup-kp(i,it)+1
            cc=0
            dcc=0
            do j=inf,isup
               phi=aper(j,1,i)+t*aper(j,2,i)
               cc=cc+per(j,i)*sin(phi)
               dcc=dcc+per(j,i)*aper(j,2,i)*cos(phi)
            enddo
            ct=cc+t*ct
            dct=dcc+t*dct
            isup=inf-1
         enddo
!
         c(i)=c(i)+ct
         dc(i)=dc(i)+dct
!
      enddo
!
      phi=al(5)
      do i=4,1,-1
         phi=phi*t+al(i)
      enddo
!
      v=c(1)+phi
      u=c(2)
      r=385000.57d0+c(3)
      v=mod(v,360.d0)
      if (v.lt.0) v=v+360.d0
!
      phi=dal(4)
!
      do i=3,1,-1
         phi=phi*t+dal(i)
      enddo
!
      vp=dc(1)+phi
      up=dc(2)
      rp=dc(3)
!
      return
      end
!
!
!
      subroutine TITRE (iprint,nul,ipr,icd,ir,idt,idf)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Impression des titres des resultats.
!
! --- Entree -----------------------------------------------------------
!
!     iprint   Indice d'impression (entier).
!              iprint = 1 pour l'Ecran.
!              iprint = 2 pour le Fichier impression.
!
!     nul      Numero d'unite logique du fichier Impression (entier).
!
!     ipr      Indice de la precision (entier).
!              ipr = 1 pour Haute precision.
!              ipr = 2 pour Moyenne precision.
!              ipr = 3 pour Faible precision.
!
!     icd      Indice des coordonnees (entier).
!              icd = 1 pour les coordonnees geometriques.
!              icd = 2 pour les coordonnees apparentes.
!              icd = 3 pour les elements orbitaux osculateurs.
!
!     ir       Indice du Repere (entier).
!              ir = 1 pour le Repere ecliptique moyen de la date
!              ir = 2 pour le Repere equatorial moyen de la date.
!              ir = 3 pour le Repere ecliptique vrai.
!              ir = 4 pour le Repere equatorial vrai.
!
!     idt      Indice du Temps (entier).
!              idt = 1 pour le Temps Dynamique.
!              idt = 2 pour le Temps Universel.
!
!     idf      Indice du type des dates (entier).
!              idt = 1 pour les Dates Juliennes.
!              idt = 2 pour les Dates Calendrier.
!
! ----------------------------------------------------------------------
!
      select case (iprint)
!
      case (1)
!
         call EFFECR
         write (*,1000)
!
         select case (icd)
         case (1)
            write (*,1001) 
!V98        write (*,1001) ipr
         case (2)
            write (*,1002) 
!V98        write (*,1002) ipr
         case (3)
            write (*,1003) 
!V98        write (*,1003) ipr
            write (*,2001)
            if (idf.eq.1) write (*,5001)
            if (idf.eq.2) write (*,5002)
            return
         endselect
!
         select case (ir)
         case (1)
            write (*,2001)
         case (2)
            write (*,2002)
         case (3)
            write (*,2003)
         case (4)
            write (*,2004)
         endselect
!
         ncas=(idf-1)*4+(idt-1)*2+2-mod(ir,2)
         select case (ncas)
         case (1)
            write (*,3001)
         case (2)
            write (*,3002)
         case (3)
            write (*,3003)
         case (4)
            write (*,3004)
         case (5)
            write (*,4001)
         case (6)
            write (*,4002)
         case (7)
            write (*,4003)
         case (8)
            write (*,4004)
         endselect
!
      case (2)
!
         select case (icd)
         case (1)
            write (nul,1001) 
!V98        write (nul,1001) ipr
         case (2)
            write (nul,1002) 
!V98        write (nul,1002) ipr
         case (3)
            write (nul,1003) 
!V98        write (nul,1003) ipr
            write (nul,2001)
            if (idf.eq.1) write (nul,5001)
            if (idf.eq.2) write (nul,5002)
            return
         endselect
!
         select case (ir)
         case (1)
            write (nul,2001)
         case (2)
            write (nul,2002)
         case (3)
            write (nul,2003)
         case (4)
            write (nul,2004)
         endselect
!
         ncas=(idf-1)*4+(idt-1)*2+2-mod(ir,2)
         select case (ncas)
         case (1)
            write (nul,3001)
         case (2)
            write (nul,3002)
         case (3)
            write (nul,3003)
         case (4)
            write (nul,3004)
         case (5)
            write (nul,4001)
         case (6)
            write (nul,4002)
         case (7)
            write (nul,4003)
         case (8)
            write (nul,4004)
         endselect
!
      endselect
!
      return
!
! --- Formats ----------------------------------------------------------
!
1000  format (2x,'JANUS GF9812',17x,'LUNE de -4000 a 8000',17x,
     .        'M&J CHAPRONT'/2x,78('-'))
!
1001  format (/2x,t30,'COORDONNEES GEOMETRIQUES')
1002  format (/2x,t31,'COORDONNEES APPARENTES')
1003  format (/2x,t27,'ELEMENTS ORBITAUX OSCULATEURS')
!V98  1001  format (/2x,t23,'COORDONNEES GEOMETRIQUES (precision',i2,')')
!V98  1002  format (/2x,t24,'COORDONNEES APPARENTES (precision',i2,')')
!V98  1003  format (/2x,t20,'ELEMENTS ORBITAUX OSCULATEURS (precision',i2,')')
!
2001  format (2x,t22,'Ecliptique et Equinoxe moyens de la date'/)
2002  format (2x,t23,'Equateur et Equinoxe moyens de la date'/)
2003  format (2x,t20,'Ecliptique moyen de la date et Equinoxe vrai'/)
2004  format (2x,t29,'Equateur et Equinoxe vrais'/)
!
3001  format (2x,t12,'Date Julienne (ET)',
     .        6x,'Longitude',5x,'Latitude',7x,'Distance')
3002  format (2x,t12,'Date Julienne (ET)',
     .        5x,'Asc droite',4x,'Declinaison',5x,'Distance')
3003  format (2x,t6,'Date Julienne (UT)',
     .        6x,'Longitude',5x,'Latitude',7x,'Distance',
     .        5x,'ET-UT')
3004  format (2x,t6,'Date Julienne (UT)',
     .        5x,'Asc droite',4x,'Declinaison',5x,'Distance',
     .        5x,'ET-UT')
!
4001  format (2x,t10,1x,'J',2x,'M',3x,'An',5x,'H (ET)',
     .        5x,'Longitude',5x,'Latitude',7x,'Distance')
4002  format (2x,t10,1x,'J',2x,'M',3x,'An',5x,'H (ET)',
     .        4x,'Asc droite',4x,'Declinaison',5x,'Distance')
4003  format (2x,t4,1x,'J',2x,'M',3x,'An',5x,'H (UT)',
     .        5x,'Longitude',5x,'Latitude',7x,'Distance',
     .        5x,'ET-UT')
4004  format (2x,t4,1x,'J',2x,'M',3x,'An',5x,'H (UT)',
     .        4x,'Asc droite',4x,'Declinaison',5x,'Distance',
     .        5x,'ET-UT')
!
5001  format (4x,'Date Julienne (ET)',
     .        8x,'a',8x,'e',5x,'gamma',
     .        3x,'L_Peri',3x,'L_Node',4x,'Lambda')
5002  format (3x,'J',2x,'M',3x,'An',5x,'H (ET)',
     .        7x,'a',8x,'e',5x,'gamma',
     .        3x,'L_Peri',3x,'L_Node',4x,'Lambda')
!
      end
!
!
!
      subroutine IMPC (d,nul,along,alat,r,ipr,dt,idt,idf)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Impression des coordonnees ecliptiques.
!
! --- Entree -----------------------------------------------------------
!
!     d        Date julienne (double reel).
!     nul      Numero d'unite logique du fichier Impression (entier).
!     ipr      Indice Precision (entier).
!     dt       Difference ET-UT en secondes (double reel).
!     idt      Indice Temps (entier).
!     idf      Indice Type de dates (entier).
!
! --- Sortie -----------------------------------------------------------
!
!     along    Longitude en degres (double reel).
!     alat     Latitude en degres (double reel).
!     r        Distance Terre-Lune en km (double reel).
!
! ----------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*22 adat
      character*50 acd
!
      id=int(dt)
!
      select case (idf)
      case (1)
         write (adat,1001) d
      case (2)
         call DJDATE (d,3,ij,im,ia,nh,nm,ns)
         write (adat,1002) ij,im,ia,nh,nm,ns
      endselect
!
      select case (idt)
      case (1)
         select case (ipr)
         case (1)
            write (acd,2001) along,alat,r
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(18:19).eq.' .') acd(18:19)='0.'
            if (acd(18:19).eq.'-.') acd(17:19)='-0.'
         case (2)
            rp=r/10.d0
            write (acd,2002) along,alat,rp
            if (acd(6:7).eq.' .') acd(6:7)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(37:37)='0'
         case (3)
            rp=r/100.d0
            write (acd,2003) along,alat,rp
            if (acd(7:8).eq.' .') acd(7:8)='0.'
            if (acd(20:21).eq.' .') acd(20:21)='0.'
            if (acd(20:21).eq.'-.') acd(19:21)='-0.'
            acd(36:37)='00'
         endselect
         write (*,2004) adat,acd
         write (nul,2004) adat,acd
      case (2)
         select case (ipr)
         case (1)
            write (acd,3001) along,alat,r,id
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(18:19).eq.' .') acd(18:19)='0.'
            if (acd(18:19).eq.'-.') acd(17:19)='-0.'
         case (2)
            rp=r/10.d0
            write (acd,3002) along,alat,rp,id
            if (acd(6:7).eq.' .') acd(6:7)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(37:37)='0'
         case (3)
            rp=r/100.d0
            write (acd,3003) along,alat,rp,id
            if (acd(7:8).eq.' .') acd(7:8)='0.'
            if (acd(20:21).eq.' .') acd(20:21)='0.'
            if (acd(20:21).eq.'-.') acd(19:21)='-0.'
            acd(36:37)='00'
         endselect
         write (*,3004) adat,acd
         write (nul,3004) adat,acd
      endselect
!
      return
!
! --- Formats ----------------------------------------------------------
!
1001  format (4x,f13.5,5x)
1002  format (i2.2,1x,i2.2,1x,i5,2x,i2.2,'h',i2.2,'m',i2.2,'s')
!
2001  format (f10.4,'d',2x,f10.4,'d',2x,f11.1,' km')
2002  format (f10.3,'d',2x,f10.3,'d',2x,f11.0,' km')
2003  format (f10.2,'d',2x,f10.2,'d',2x,f10.0,'  km')
2004  format (t10,a22,2x,a40)
!
3001  format (f10.4,'d',2x,f10.4,'d',2x,f11.1,' km',2x,i6,' s')
3002  format (f10.3,'d',2x,f10.3,'d',2x,f11.0,' km',2x,i6,' s')
3003  format (f10.2,'d',2x,f10.2,'d',2x,f10.0,'  km',2x,i6,' s')
3004  format (t4,a22,2x,a50)
!
      end
!
!
!
      subroutine IMPQ (d,nul,ra,de,r,ipr,dt,idt,idf)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Impression des coordonnees equatoriales
!
! --- Entree -----------------------------------------------------------
!
!     d        Date julienne (double reel).
!     nul      Numero d'unite logique du fichier Impression (entier).
!     ipr      Indice Precision (entier).
!     dt       Difference ET-UT en secondes (double reel).
!     idt      Indice Temps (entier).
!     idf      Indice Type de dates (entier).
!
! --- Sortie -----------------------------------------------------------
!
!     ra       Acension droite en heures (double reel).
!     de       Declinaison en degres (double reel).
!     r        Distance Terre-Lune en km (double reel).
!
! ----------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*22 adat
      character*50 acd
!
      id=int(dt)
!
      select case (idf)
      case (1)
         write (adat,1001) d
      case (2)
         call DJDATE (d,3,ij,im,ia,nh,nm,ns)
         write (adat,1002) ij,im,ia,nh,nm,ns
      endselect
!
      select case (idt)
      case (1)
         select case (ipr)
         case (1)
            write (acd,2001) ra,de,r
            if (acd(4:5).eq.' .') acd(4:5)='0.'
            if (acd(18:19).eq.' .') acd(18:19)='0.'
            if (acd(18:19).eq.'-.') acd(17:19)='-0.'
         case (2)
            rp=r/10.d0
            write (acd,2002) ra,de,rp
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(37:37)='0'
         case (3)
            rp=r/100.d0
            write (acd,2003) ra,de,rp
            if (acd(6:7).eq.' .') acd(6:7)='0.'
            if (acd(20:21).eq.' .') acd(20:21)='0.'
            if (acd(20:21).eq.'-.') acd(19:21)='-0.'
            acd(36:37)='00'
         endselect
         write (*,2004) adat,acd
         write (nul,2004) adat,acd
      case (2)
         select case (ipr)
         case (1)
            write (acd,3001) ra,de,r,id
            if (acd(4:5).eq.' .') acd(4:5)='0.'
            if (acd(18:19).eq.' .') acd(18:19)='0.'
            if (acd(18:19).eq.'-.') acd(17:19)='-0.'
         case (2)
            rp=r/10.d0
            write (acd,3002) ra,de,rp,id
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(37:37)='0'
         case (3)
            rp=r/100.d0
            write (acd,3003) ra,de,rp,id
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(36:37)='00'
         endselect
         write (*,3004) adat,acd
         write (nul,3004) adat,acd
      endselect
!
      return
!
! --- Formats ----------------------------------------------------------
!
1001  format (4x,f13.5,5x)
1002  format (i2.2,1x,i2.2,1x,i5,2x,i2.2,'h',i2.2,'m',i2.2,'s')
!
2001  format (f10.5,'h',2x,f10.4,'d',2x,f11.1,' km')
2002  format (f10.4,'h',2x,f10.3,'d',2x,f11.0,' km')
2003  format (f10.3,'h',2x,f10.2,'d',2x,f10.0,'  km')
2004  format (t10,a22,2x,a40)
!
3001  format (f10.5,'h',2x,f10.4,'d',2x,f11.1,' km',2x,i6,' s')
3002  format (f10.4,'h',2x,f10.3,'d',2x,f11.0,' km',2x,i6,' s')
3003  format (f10.3,'h',2x,f10.2,'d',2x,f10.0,'  km',2x,i6,' s')
3004  format (t4,a22,2x,a50)
!
        end
!
!
!
      subroutine IMPS (d,nul,a,e,gam,pom,gom,alam,ipr,idf)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Impression des elements orbitaux osculateurs.
!
! --- Entree -----------------------------------------------------------
!
!     d        Date julienne (double reel).
!     nul      Numero d'unite logique du fichier Impression (entier).
!     ipr      Indice Precision (entier).
!     idf      Indice Type de dates (entier).
!
! --- Sortie -----------------------------------------------------------
!
!     a        Demi-grand axe en km (double reel).
!     e        Exentricite (double reel).
!     gam      Sinus de la demi inclinaison (double reel).
!     pom      Longitude du perigee en degres (double reel).
!     gom      Longitude du noeud en degres (double reel).
!     alam     Longitude moyenne en degres (double reel).
!
! ----------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*22 adat
      character*54 acd
!
      select case (idf)
      case (1)
         write (adat,1001) d
      case (2)
         call DJDATE (d,3,ij,im,ia,nh,nm,ns)
         write (adat,1002) ij,im,ia,nh,nm,ns
      endselect
!
      select case (ipr)
      case (1)
         write (acd,2001) a,e,gam,pom,gom,alam
         acd(7:7)='k'
         acd(10:10)='0'
         acd(18:18)='0'
         if (acd(30:31).eq.' .') acd(30:31)='0.'
         if (acd(38:39).eq.' .') acd(38:39)='0.'
         if (acd(47:48).eq.' .') acd(47:48)='0.'
      case (2)
         ap=a/100.d0
         write (acd,2002) ap,e,gam,pom,gom,alam
         acd(5:7)='00k'
         acd(11:11)='0'
         acd(19:19)='0'
         if (acd(31:32).eq.' .') acd(31:32)='0.'
         if (acd(39:40).eq.' .') acd(39:40)='0.'
         if (acd(49:50).eq.' .') acd(49:50)='0.'
      case (3)
         ap=a/100.d0
         write (acd,2003) ap,e,gam,pom,gom,alam
         acd(5:7)='00k'
         acd(12:12)='0'
         acd(20:20)='0'
         acd(34:34)='d'
         if (acd(40:41).eq.' .') acd(40:41)='0.'
         if (acd(50:51).eq.' .') acd(50:51)='0.'
      endselect
!
      write (*,3000) adat,acd
      write (nul,3000) adat,acd
!
      return
!
! --- Formats ----------------------------------------------------------
!
1001  format (4x,f13.5,5x)
1002  format (i2.2,1x,i2.2,1x,i5,2x,i2.2,'h',i2.2,'m',i2.2,'s')
!
2001  format (f7.0,'m',1x,f7.5,1x,f8.6,2x,
     .        f6.2,'d',1x,f7.3,'d',1x,f8.4,'d')
2002  format (f5.0,'  m',1x,f7.4,1x,f8.5,2x,
     .        f6.1,'d',1x,f7.2,'d',1x,f8.2,'d')
2003  format (f5.0,'  m',1x,f7.3,1x,f8.4,2x,
     .        f7.0,1x,f7.1,'d',1x,f8.2,'d')
!
3000  format (2x,a22,2x,a54)
!
        end
!
!
!
      subroutine DCDJ (cd,d,ierr,datcal)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Changement d'une Date Calendrier en Date Julienne.
!
! --- Entree -----------------------------------------------------------
!
!     cd       Date Calendrier +/-aaaammjj.hhmmss (double reel).
!
! --- Sortie -----------------------------------------------------------
!
!     d        Date Julienne (double reel).
!
!     ierr     Indice Erreur (entier).
!              ierr = 0  pas d'erreur.
!              ierr = 1  date entre le 05/10/1582 et la 15/10/1582.
!              ierr = 2  date non correcte.
!
!     datcal   Date Calendrier si ierr=0 ou 1 (character*30).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*4 cmo(12)
      character*16 car
      character*30 datcal
!
! --- Initialisation ---------------------------------------------------
!
      data cmo/
     .'Janv','Fevr','Mars','Avr.','Mai','Juin',
     .'Juil','Aout','Sept','Oct.','Nov.','Dec.'/
!
      d=0.d0
      ierr=0
      datcal=' '
!
! --- Decodage de la Date Calendrier -----------------------------------
!
      write (car,1000) cd
      if (car(10:10).ne.'.') then
         ierr=2
         return
      endif
!
      read (car(1:5),1001) iy
      read (car(6:7),1002) im
      read (car(8:9),1002) id
      read (car(11:12),1002) ih
      read (car(13:14),1002) imn
      read (car(15:16),1002) is
      read (car(1:9),1003) iw
!
! --- Verification des parametres --------------------------------------
!
      if (im.lt.1.or.im.gt.12.or.
     .    id.gt.31.or.ih.gt.24.or.imn.gt.59.or.is.gt.59) then
          ierr=2
          return
      endif
!
      write (datcal,2000) id,cmo(im),iy,ih,imn,is
!
! --- Calcul de la Date Julienne ---------------------------------------
!
      if (im.le.2) then
         ims=im+13
         iyp=iy-1
      else
         ims=im+1
         iyp=iy
      endif
!
      if (iw.ge.15821015) then
         ia=iyp/100
         ib=2-ia+ia/4
      else
         if (iw.ge.15821005) then
            ierr=1
            return
         endif
         ib=0
      endif
!
      w=365.25d0*iyp
      if (iyp.lt.0) w=w-0.75d0
!
      d=int(w)+int(30.60001d0*ims)+id+ib+ih/24.d0+imn/1440.d0
     . +is/86400.d0+1720994.5d0
!
      return
!
! --- Formats ----------------------------------------------------------
!
1000  format(f16.6)
1001  format(i5)
1002  format(i2)
1003  format(i9)
2000  format(i2.2,1x,a4,1x,i5,3x,i2.2,'h',i2.2,'m',i2.2,'s')
!
        end
!
!
!
      subroutine DJDATE (tdj,icode,jour,mois,kan,nh,nm,ns)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Changement d'une Date Julienne en Date Calendrier.
!
! --- Entree ------------------------------------------------------------
!
!     tdj         Date julienne (reel dp).
!
!     icode       Indice code de la conversion (entier)
!                 icode=0 conversion tdj -> jour-mois-kan et 00-00-00.
!                 icode=1 conversion tdj -> jour-mois-kan et nh-00-00.
!                 icode=2 conversion tdj -> jour-mois-kan et nh-nm-00.
!                 icode=3 conversion tdj -> jour-mois-kan et nh-nm-ns.
!
! --- Sortie -----------------------------------------------------------
!
!     jour        Quantieme du mois (entier).
!     mois        Numero du mois (entier).
!     kan         Millesime astronomique de l'annee (entier).
!
!     nh          Nombre d'heures (entier).
!     nm          Nombre de minutes (entier).
!     ns          Nombre de secondes (entier).
!
! --- Remarque ---------------------------------------------------------
!
!     Date du calendier (Gregorien/Julien) :
!     Le calendier Gregorien commence le 15/10/1582 a 0h (DJ2299160.5).
!     Avant cette date on utilise le calendrier Julien.
!
!     On utilise le formulaire de J. Meeus.
!     La date julienne tdj doit etre superieure ou egale a 0.
!     Si tdj<0, alors jour=mois=kan=nh=nm=ns=0.
!
!     L'annee kan est comptee en millesimes astronomiques.
!     L'annee 1 avant J.C. des historiens correspond a l'an 0 des
!     astronomes et ainsi de suite.
!     Ex : l'an 80 avant J.C. correspond a l'an astronomique -79.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
! --- Initialisation ---------------------------------------------------
!
      jour=0
      mois=0
      kan=0
      nh=0
      nm=0
      ns=0
!
      if (tdj.lt.0.d0.or.icode.lt.0.or.icode.gt.3) return
      if (icode.eq.0) dt=0.0d0
      if (icode.eq.1) dt=0.5d0/24.d0
      if (icode.eq.2) dt=0.5d0/1440.d0
      if (icode.eq.3) dt=0.5d0/86400.d0
!
! --- Conversion -------------------------------------------------------
!
      t=tdj+dt+0.5d0
      z=int(t)
      f=t-z
!
      if (z.lt.2299161.d0) then
         a=z
      else
         alpha=int((z-1867216.25d0)/36524.25d0)
         a=z+1.d0+alpha-int(alpha/4.d0)
      endif
!
      b=a+1524.d0
      c=int((b-122.1d0)/365.25d0)
      d=int(365.25d0*c)
      e=int((b-d)/30.6001d0)
!
      jour=b-d-int(30.6001d0*e)
!
      if (e.lt.13.5d0) then
         mois=e-1.d0
      else
         mois=e-13.d0
      endif
!
      if (mois.lt.2.5d0) then
         kan=c-4715.d0
      else
         kan=c-4716.d0
      endif
!
      f=f*24.d0
      nh=f
      f=(f-nh)*60.d0
      nm=f
      f=(f-nm)*60.d0
      ns=f
!
! --- Sortie -----------------------------------------------------------
!
      return
      end
!
!
!
      subroutine UT (dj,dt,ierr,limut,tlimut)
! ----------------------------------------------------------------------
!
!     Ref : GF-0212
!
! --- Objet ------------------------------------------------------------
!
!     Difference ET-UT.
!
! --- Entree -----------------------------------------------------------
!
!     dj       Date Julienne (double reel).
!
! --- Sortie -----------------------------------------------------------
!
!     dt       ET-UT en secondes (double reel).
!
!     ierr     Indice Erreur (entier).
!              ierr = 0  Pas d'erreur.
!              ierr = 1  Erreur : Date posterieure a la limite UT.
!
!     limut    Annee limite UT (entier).
!     tlimut   Date Julienne limite UT (entier).
!
! --- Remarques --------------------------------------------------------
!
!     La date dj peut etre comptee indifferemment en ET ou en UT.
!
!     Pour les periodes passees, on utilise les formules de Stepheson
!     et Morrison (1984) en tenant compte de l'acceleration seculaire 
!     liee aux termes de marees de la solution ELP2000-85.

!     Pour la periode contemporaine, on utilise les valeurs observees
!     de TAI-UT1. Les valeurs de ET-TU sont exprimees en secondes dans 
!     le tableau idtsm. Elles sont donnees tous les 5 ans entre 1600 et 
!     1780 et chaque annee a partir de 1780 (voir mise a jour annuelle).

!     Au dela de l'annnee limite de mise a jour il n'y a pas de calcul
!     de ET-UT (dt=0.0).
!
! --- Mise a jour annuelle ---------------------------------------------
!
!     - augmenter de 1 la taille du tableau idtsm (parameter idim),
!     - ajouter la valeur de ET-TU en secondes au tableau idtsm.
!
! ----------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
!***********************************************************************
!     Derniere mise a jour : jusqu'a 2005,  ET-TU=64 secondes (idim=262)
      parameter (idim=262)
      dimension idtsm(idim)
      data idtsm/
     .  80,78,75,73,70,68,65,52,43,36,30,26,20,16,10,
     .  06,01,-2,-3,-4,-3,-3,-2,-1,00,01,01,03,03,05,
     .  05,06,08,09,09,11,11,
     .  11,11,11,11,12,12,12,12,12,12,12,11,11,11,11,
     .  10,10,10,09,09,09,09,08,08,08,08,08,08,08,08,
     .  09,09,09,09,09,09,09,09,09,09,08,08,08,07,07,
     .  06,06,06,05,05,04,04,03,03,03,03,03,03,03,03,
     .  03,03,04,04,04,04,04,05,05,05,05,05,05,06,06,
     .  06,06,06,06,06,06,06,05,05,04,04,03,01,00,00,
     .  -1,-2,-3,-4,-4,-5,-6,-6,-6,-6,-6,-6,-6,-6,-7,
     .  -7,-7,-7,-6,-7,-7,-7,-7,-7,-7,-7,-6,-5,-4,-3,
     .  -2,-1,01,02,03,05,06,07,09,10,11,13,14,16,17,
     .  18,19,20,21,21,22,22,23,23,23,24,24,24,24,24,
     .  24,24,24,24,24,24,24,24,24,24,25,25,26,26,27,
     .  27,28,28,29,29,30,30,30,31,31,31,32,32,33,33,
     .  34,34,34,35,36,37,37,38,39,40,41,42,43,44,45,
     .  46,48,49,50,51,51,52,53,54,54,55,55,56,56,57,
     .  58,58,59,60,61,62,63,63,64,64,64,64,64,64,64/
!***********************************************************************
!
      dimension f1(3),f2(3)
      data f1/2177.d0,495.d0,42.4d0/
      data f2/102.d0,100.d0,23.6d0/
!
      data d2000/2451545.d0/
      data d1780/2371191.5d0/
      data sjul/36525.d0/
!
      data t948/-10.519658d0/
      data t1600/-3.999932d0/
      data t1780/-2.199959d0/
!
      ierr=0
      dt=0.d0
!
      ia=idim-36
      limut=1780+ia
      tlimut=d1780+ia*sjul/100.d0
!
      if (dj.ge.tlimut) then
         ierr=1
         return
      endif
!
      t=(dj-d2000)/sjul
!
      if ((t-t948).lt.0) then
         dt=(f1(3)*t+f1(2))*t+f1(1)
      else
         if ((t-t1600).lt.0) then
            dt=(f2(3)*t+f2(2))*t+f2(1)
         else
            if ((t-t1780).lt.0) then
               i=int((t-t1600)*20)+1
            else
               i=int((t-t1780)*100)+37
            endif
            dt=idtsm(i)
         endif
      endif
!
      return
      end
!
!
!
      subroutine EFFECR
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Objet ------------------------------------------------------------
!
!     Effacement de l'ecran et positionnement curseur en haut de page.
!
!     Selon le systeme, on utilise la commande 'cls' ou 'clear'.
!
! ----------------------------------------------------------------------
!
      call SYSTEM ('cls')
!!!   call SYSTEM ('clear')
      write (*,'(/)')
!
      return
      end
      